﻿using System;
/// <summary>
/// 矩阵
/// </summary>
public class Mat {
    /*--------------------------------------class comment
     Version   :  V1.0
     Coded by  :  syz
     Date      :  2011-06-16 16:53:40          *星期四*
    ----------------------------------------------------
    Desciption :
     *       矩阵类
     *       该类中的部分算符重载需要用到VEC类
    parameters :
     * 
    Methods    :
     *    Mat------构造函数
     *    +  --------两矩阵相加（对应元素分别相加）
     *    -  --------矩阵相减（对应元素分别相减）
     *    *  --------矩阵元素分别相乘（相当于MATLAB中的  .*算符）
     *    /  --------矩阵元素分别相除（相当与MATLAB中的  ./算符）
     * 
     *    +  --------矩阵加实数（所有元素分别加该实数）
     *    -  --------矩阵减实数（所有元素分别减该实数）
     *    *  --------矩阵乘实数（所有元素分别乘该实数）
     *    /  --------矩阵除实数（所有元素分别除该实数）
     *    
     *               以下再次重载四则算符，
     *               如此可以在矩阵左右均可进行与实数的双目运算
     *    +  --------实数加矩阵（所有元素分别加该实数）
     *    -  --------实数减矩阵（所有元素分别减该实数）
     *    *  --------实数乘矩阵（所有元素分别乘该实数）
     *    /  --------<无>---不重载--(即不能用实数除以矩阵)
     *    
     *    |  --------数学意义下的矩阵相乘
     *    |  --------矩阵与向量相乘（线性变换）
     *    ~  --------矩阵转置（向量中~该算符用于求2范数平方）
    --------------------------------------------------*/

    public int dim1,dim2;  //数组维数
    public double[,] ele;

    public Mat(int m,int n) {
        //构造函数
        dim1 = m;
        dim2 = n;
        //矩阵维数

        ele = new double[m,n];
        //用二维数组构造数学意义下的矩阵
        //矩阵元素保存于对象ele中
    }


    //-------两个矩阵加法算符重载
    //------矩阵元素分别相加
    public static Mat operator +(Mat a1,Mat a2) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a3 = new Mat(m,n);
        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a3.ele[i,j] = a1.ele[i,j] + a2.ele[i,j];

        return a3;

    }


    //-------两个矩阵减法算符重载
    //------矩阵元素分别相减
    public static Mat operator -(Mat a1,Mat a2) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a3 = new Mat(m,n);
        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a3.ele[i,j] = a1.ele[i,j] - a2.ele[i,j];

        return a3;

    }

    //-------两个矩阵乘法算符重载
    //------矩阵元素分别相乘，相当于MATLAB中的   .*
    // 要求两个矩阵维数相同，矩阵类不进行个数判断
    public static Mat operator *(Mat a1,Mat a2) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a3 = new Mat(m,n);
        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a3.ele[i,j] = a1.ele[i,j] * a2.ele[i,j];

        return a3;

    }

    //-------两个矩阵除法算符重载
    //------矩阵元素分别相除，相当于MATLAB中的   ./
    // 要求两个矩阵维数相同，矩阵类不进行个数判断
    public static Mat operator /(Mat a1,Mat a2) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a3 = new Mat(m,n);
        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a3.ele[i,j] = a1.ele[i,j] / a2.ele[i,j];

        return a3;
    }

    //矩阵加实数算符重载
    //各分量分别加实数
    public static Mat operator +(Mat a1,double x) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a2 = new Mat(m,n);

        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a2.ele[i,j] = a1.ele[i,j] + x;

        return a2;

    }


    //矩阵减实数算符重载
    //各分量分别减实数
    public static Mat operator -(Mat a1,double x) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a2 = new Mat(m,n);

        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a2.ele[i,j] = a1.ele[i,j] - x;

        return a2;
    }

    //矩阵乘以实数算符重载
    //各分量分别乘以实数
    public static Mat operator *(Mat a1,double x) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a2 = new Mat(m,n);

        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a2.ele[i,j] = a1.ele[i,j] * x;

        return a2;
    }


    //矩阵除以实数算符重载
    //各分量分别除以实数
    public static Mat operator /(Mat a1,double x) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a2 = new Mat(m,n);

        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a2.ele[i,j] = a1.ele[i,j] / x;

        return a2;
    }


    //实数加矩阵算符重载
    //各分量分别加实数
    public static Mat operator +(double x,Mat a1) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a2 = new Mat(m,n);

        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a2.ele[i,j] = a1.ele[i,j] + x;

        return a2;

    }


    //实数减矩阵算符重载
    //各分量分别减实数
    public static Mat operator -(double x,Mat a1) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a2 = new Mat(m,n);

        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a2.ele[i,j] = a1.ele[i,j] - x;

        return a2;
    }

    //实数乘矩阵算符重载
    //各分量分别乘以实数
    public static Mat operator *(double x,Mat a1) {
        int m,n;
        m = a1.dim1;
        n = a1.dim2;

        Mat a2 = new Mat(m,n);

        for (int i = 0 ; i < m ; i++)
            for (int j = 0 ; j < n ; j++)
                a2.ele[i,j] = a1.ele[i,j] * x;

        return a2;
    }



    //数学上的矩阵相乘
    public static Mat operator |(Mat a1,Mat a2) {
        int m,n,p,q;

        m = a1.dim1;
        n = a1.dim2;

        p = a2.dim1;
        q = a2.dim2;

        if (n != p) throw new Exception("Inner matrix dimensions must agree！");
        //如果矩阵维数不匹配给出告警信息

        //新矩阵，用于存放结果
        Mat a3 = new Mat(m,q);

        int i,j;


        for (i = 0 ; i < m ; i++)
            for (j = 0 ; j < q ; j++) {
                a3.ele[i,j] = 0.0;
                for (int k = 0 ; k < n ; k++)
                    a3.ele[i,j] = a3.ele[i,j] + a1.ele[i,k] * a2.ele[k,j];
            }

        return a3;
    }

    //矩阵乘以向量(线性变换）
    //即 b=Ax
    public static Vec operator |(Mat A,Vec x) {
        int m,n,p;
        m = A.dim1;
        n = A.dim2;

        p = x.dim;

        if (n != p) throw new Exception("Inner matrix dimensions must agree！");
        //如果矩阵维数不匹配，给出告警信息

        Vec b = new Vec(m);

        for (int i = 0 ; i < m ; i++) {
            b.ele[i] = 0.0;

            for (int k = 0 ; k < n ; k++)
                b.ele[i] = b.ele[i] + A.ele[i,k] * x.ele[k];
        }

        return b;
    }


    //矩阵转置
    public static Mat operator ~(Mat A) {
        int m,n;
        m = A.dim1;
        n = A.dim2;

        Mat TA = new Mat(n,m);
        for (int i = 0 ; i < n ; i++)
            for (int j = 0 ; j < m ; j++)
                TA.ele[i,j] = A.ele[j,i];

        return TA;

    }


    public double this[int index1,int index2]
        // 矩阵索引器
    {
        get {
            return ele[index1,index2];
        }
        set {
            ele[index1,index2] = value;
        }
    }


}

/// <summary>
/// 向量
/// </summary>
public class Vec
    /*--------------------------------------class comment
     Version   :  V1.0
     Coded by  :  syz
     Date      :  2011-06-15 12:50:11          *星期三*
    ----------------------------------------------------
    Desciption :
                向量类
    parameters :
     *   dim------向量维数
     *   ele------向量元素
     * 
    Methods    :
     *    =    -----------向量相等
     *                    注意该方法不需要自行编写重载函数，
     *                    c#已经实现了对类的等号算符重载，并禁止
     *                    用户编写等号算符（=）重载函数
     *    +    -----------两向量对应元素相加
     *    -    -----------两向量对应元素相减
     *    *    -----------两向量对应元素相乘（相当于MATLAB中的  .*算符）
     *    /    -----------两向量对应元素相除（相当于MATLAB中的  ./算符）
     *    
     *    +  -----------重载向量加一实数（所有元素）
     *    -  -----------重载向量减一实数（所有元素）
     *    *  -----------重载向量乘一实数（所有元素）
     *    /  -----------重载向量除一实数（所有元素）
     *   
     *    +  -----------重载实数加向量（所有元素）
     *    -  -----------重载实数减向量（所有元素）
     *    *  -----------重载实数乘以向量（所有元素）
     *    /  -----------<无>--不重载--<即不能用实数除以向量>
     * 
     *    |  （VEC1|ve2）----向量内积
     *    ~  （~VEC1）-------向量2范数的平方
     *    -  （-VEC1）-------向量取负
     *    ^  （VEC1^VEC2）---向量外积
     * 
    --------------------------------------------------*/
{
    public int dim;  //数组维数
    public double[] ele;

    public Vec(int m) {
        //构造函数
        dim = m;
        ele = new double[dim];
        //用一维数组构造向量
    }


    //-------两个向量加法算符重载
    //------分量分别相加
    public static Vec operator +(Vec v1,Vec v2) {
        int N0;

        //获取变量维数
        N0 = v1.dim;

        Vec v3 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v3.ele[j] = v1.ele[j] + v2.ele[j];
        }
        return v3;
    }

    //-----------------向量减法算符重载
    //-----------------分量分别想减
    public static Vec operator -(Vec v1,Vec v2) {
        int N0;
        //获取变量维数
        N0 = v1.dim;

        Vec v3 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v3.ele[j] = v1.ele[j] - v2.ele[j];
        }
        return v3;
    }

    //----------------向量乘法算符重载
    //-------分量分别相乘，相当于MATLAB中的  .*算符
    public static Vec operator *(Vec v1,Vec v2) {

        int N0;

        //获取变量维数
        N0 = v1.dim;

        Vec v3 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v3.ele[j] = v1.ele[j] * v2.ele[j];
        }
        return v3;
    }


    //---------------向量除法算符重载
    //--------------分量分别相除，相当于MATLAB中的   ./算符
    public static Vec operator /(Vec v1,Vec v2) {
        int N0;

        //获取变量维数
        N0 = v1.dim;

        Vec v3 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v3.ele[j] = v1.ele[j] / v2.ele[j];
        }
        return v3;
    }


    //向量减加实数
    //各分量分别加实数
    public static Vec operator +(Vec v1,double a) {
        //向量数加算符重载
        int N0;

        //获取变量维数
        N0 = v1.dim;

        Vec v2 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v2.ele[j] = v1.ele[j] + a;
        }
        return v2;
    }

    //向量减实数
    //各分量分别减实数
    public static Vec operator -(Vec v1,double a) {
        //向量数加算符重载
        int N0;

        //获取变量维数
        N0 = v1.dim;

        Vec v2 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v2.ele[j] = v1.ele[j] - a;
        }
        return v2;
    }


    //向量 数乘
    //各分量分别乘以实数
    public static Vec operator *(Vec v1,double a) {
        int N0;

        //获取变量维数
        N0 = v1.dim;

        Vec v2 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v2.ele[j] = v1.ele[j] * a;
        }
        return v2;
    }


    //向量 数除
    //各分量分别除以实数
    public static Vec operator /(Vec v1,double a) {
        int N0;

        //获取变量维数
        N0 = v1.dim;

        Vec v2 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v2.ele[j] = v1.ele[j] / a;
        }
        return v2;
    }


    //实数加向量
    public static Vec operator +(double a,Vec v1) {
        //向量数加算符重载
        int N0;

        //获取变量维数
        N0 = v1.dim;

        Vec v2 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v2.ele[j] = v1.ele[j] + a;
        }
        return v2;
    }

    //实数减向量
    public static Vec operator -(double a,Vec v1) {
        //向量数加算符重载
        int N0;

        //获取变量维数
        N0 = v1.dim;

        Vec v2 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v2.ele[j] = v1.ele[j] - a;
        }
        return v2;
    }


    //向量 数乘
    public static Vec operator *(double a,Vec v1) {
        int N0;

        //获取变量维数
        N0 = v1.dim;

        Vec v2 = new Vec(N0);

        int j;
        for (j = 0 ; j < N0 ; j++) {
            v2.ele[j] = v1.ele[j] * a;
        }
        return v2;
    }



    //---------------向量内积
    public static double operator |(Vec v1,Vec v2) {
        int N0,M0;

        //获取变量维数
        N0 = v1.dim;

        M0 = v2.dim;

        if (N0 != M0)
            throw new Exception("Inner vector dimensions must agree！");
        //如果向量维数不匹配，给出告警信息

        double sum;
        sum = 0.0;

        int j;
        for (j = 0 ; j < N0 ; j++) {
            sum = sum + v1.ele[j] * v2.ele[j];
        }
        return sum;
    }


    //-----------向量外积（相当于列向量，乘以横向量）
    public static Mat operator ^(Vec v1,Vec v2) {
        int N0,M0;

        //获取变量维数
        N0 = v1.dim;

        M0 = v2.dim;

        if (N0 != M0)
            throw new Exception("Inner vector dimensions must agree！");
        //如果向量维数不匹配，给出告警信息

        Mat vvmat = new Mat(N0,N0);

        for (int i = 0 ; i < N0 ; i++)
            for (int j = 0 ; j < N0 ; j++)
                vvmat[i,j] = v1[i] * v2[j];

        //返回外积矩阵
        return vvmat;

    }


    //---------------向量模的平方
    public static double operator ~(Vec v1) {
        int N0;

        //获取变量维数
        N0 = v1.dim;

        double sum;
        sum = 0.0;

        int j;
        for (j = 0 ; j < N0 ; j++) {
            sum = sum + v1.ele[j] * v1.ele[j];
        }
        return sum;
    }

    //  负向量 
    public static Vec operator -(Vec v1) {
        int N0 = v1.dim;

        Vec v2 = new Vec(N0);

        for (int i = 0 ; i < N0 ; i++)
            v2.ele[i] = -v1.ele[i];

        return v2;
    }

    public double this[int index]
        /*------------------------------------------ comment
        Author      : syz 
        Date        : 2011-07-03 19:53:28         
        ---------------------------------------------------
        Desciption  : 创建索引器
         *
        Post Script :
         *
        paramtegers :
         *
        -------------------------------------------------*/
    {
        get {
            //get accessor
            return ele[index];

        }
        set {
            //set accessor
            ele[index] = value;

        }

    }
}
